function fsb_diag(idat,imean,sandbox,inparg)

% FSB : Calculate and display diagnostic maps
%
% EXAMPLE:
% fsb_diag(maptype,idat,imean,sandbox,slice_x,slice_y,slicez)
%
% INPUT:
% inparg.maptype: type of map selected
% idat: 4-D image data
% imean: mean image data
% sandbox:      sandbox experiment struct
% inparg.slice_n: Selected slices
% inparg.hrf_pred:     Selected hemodynamic predictor
%
% OUTPUT:
% Diagnostic maps
%
% CALLED BY:
%
% FSB.m
% fsb_functional_map.m
%
% NOTES:
%
% Copyright 2010 MPI for Biological Cybernetics
% Author: Steffen Stoewer
% License:GNU GPL, no express or implied warranties
%
%$ Revision 1.0
%
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

[a b c d] = size(idat);

switch inparg.maptype;

    case 'Scaninfo'
        inparg.hnoise = 6;
        inparg.maptypenum = 1;
        Disp_scaninfo
    case 'Plot slices'
        inparg.hnoise = 6;
        inparg.maptypenum = 2;
        Disp_slices
    case 'Voxel Time Course';
        inparg.hnoise = 6;
        inparg.maptypenum = 3;
        fsb_plot_voxel_timecourse(sandbox,idat,inparg.slice_n);
    case 'SNR Map'
        inparg.maptypenum = 4;
        if inparg.hnoise>5
            inparg.hnoise = 4;
            inparg.noiseplot = [size(idat,2)-3 3];
            inparg.noise_origin = [3 3];
            inparg.noise = [3 3 3 size(idat,2)-3];
            inparg.noisetd = [size(idat,3)-3 size(idat,3)];

        end
        Diag_SNR;
    case 'ROI SNR Map'
        inparg.maptypenum = 5;
        inparg.hnoise = 5;
        inparg.noise_origin(1) = inparg.slice_n(2)-inparg.roi;
        inparg.noise_origin(2) = inparg.slice_n(1)-inparg.roi;
        inparg.noiseplot(1) = inparg.roi*2;
        inparg.noiseplot(2) = inparg.roi*2;
        inparg.noise = [inparg.slice_n(2)-inparg.roi inparg.slice_n(2)+inparg.roi inparg.slice_n(1)-inparg.roi inparg.slice_n(1)+inparg.roi];
        inparg.noisetd = [inparg.slice_n(3)-inparg.roi inparg.slice_n(3)+inparg.roi];
        Diag_SNR;
    case 'Sandmap'
        inparg.hnoise = 6;
        inparg.maptypenum = 14;
        Diag_Sandmap;
    case 'STD Map'
        inparg.hnoise = 6;
        inparg.maptypenum = 6;
        Diag_STD
    case 'STD Map2'
        inparg.hnoise = 6;
        inparg.maptypenum = 13;
        Diag_STD2
    case 'TN Map'
        inparg.hnoise = 6;
        inparg.maptypenum = 7;
        Diag_TN
    case 'CV Map'
        inparg.hnoise = 6;
        inparg.maptypenum = 8;
        Diag_CV
    case 'Z-score Map';
        inparg.hnoise = 6;
        inparg.maptypenum = 9;
        Diag_zscore
    case 'Special Map';
        inparg.hnoise = 6;
        inparg.maptypenum = 10;
        Diag_special
    case 'AB Map';
        inparg.hnoise = 6;
        inparg.maptypenum = 11;
        Diag_AB
    case 'Ghost Map';
        inparg.hnoise = 6;
        inparg.maptypenum = 12;
        Diag_Ghost
    case 'Noise Time Course';
        inparg.hnoise = 6;
        inparg.maptypenum = 12;
        fsb_noise_map(sandbox,inparg.hrf_pred);
    case 'Volume Histogram';
        inparg.hnoise = 6;
        inparg.maptypenum = 14;
        Diag_Hist
    case 'Euclidian Distance';
        inparg.hnoise = 6;
        inparg.maptypenum = 15;
        fsb_diag_euclid(idat,1);

end


%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Show brain slices only
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    function Disp_slices

        idat_map = single(idat(:,:,:,inparg.slice_n(4)));
        inparg.map_txt = 'Brain slices';
        fsb_functional_map(idat_map,inparg,sandbox);

    end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Display scan info
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    function Disp_scaninfo
        
        

        try
            fm_par = FM_calc_bw(sandbox);
            disp(' ')
            disp('=================================');
            disp(['Scan information for scan ' num2str(sandbox.scanpar.acqp.expno)]);
            disp('=================================');
            disp(['Protocol name: ' sandbox.scanpar.acqp.ACQ_protocol_name]);
            disp(['Pulse Program: ' sandbox.scanpar.acqp.PULPROG]);
            disp(['Acquisition Method: ' sandbox.scanpar.acqp.ACQ_method]);
            disp(['Matrix size: ' num2str(sandbox.scanpar.acqp.PVM_Matrix(1)) 'x' num2str(sandbox.scanpar.acqp.PVM_Matrix(2))]);
            disp(['Acquisition resolution: ' num2str(sandbox.scanpar.acqp.PVM_SpatResol(1)) 'x' num2str(sandbox.scanpar.acqp.PVM_SpatResol(2))]);
            disp(['Number of slices: ' num2str(sandbox.scanpar.acqp.NSLICES)]);
            disp(['Slice offset: ' num2str(sandbox.scanpar.acqp.ACQ_slice_offset(1)) ' ' num2str(sandbox.scanpar.acqp.ACQ_slice_offset(end))]);

            if sandbox.scanpar.acqp.ACQ_slice_thick<10
                disp(['Slice thickness: ' num2str(sandbox.scanpar.acqp.ACQ_slice_thick)]);
                disp(['Slice center: ' num2str(find(sandbox.scanpar.acqp.ACQ_slice_offset==0))]);
            else
                disp(['Slice thickness PVM: ' num2str(sandbox.scanpar.acqp.PVM_SpatResol(3))]);
                stack_height = sandbox.scanpar.acqp.ACQ_slice_thick;
                slice_center = stack_height/2+sandbox.scanpar.acqp.ACQ_slice_offset;
                disp(['Slice center: ' num2str(slice_center)]);
            end

            disp(['Reco size: ' num2str(sandbox.scanpar.reco.RECO_size(1)) 'x' num2str(sandbox.scanpar.reco.RECO_size(2)) ]);
            disp(['Reco FOV: ' num2str(sandbox.scanpar.reco.RECO_fov(1)) 'x' num2str(sandbox.scanpar.reco.RECO_fov(2))]);
            disp(['Echo Time (TE): ' num2str(sandbox.scanpar.acqp.ACQ_echo_time)]);
            disp(['Repetition Time (TR): ' num2str(sandbox.scanpar.acqp.ACQ_repetition_time)]);
            disp(' ')
            fm_par = FM_calc_bw(sandbox);
        catch
            disp('Full scan information not available')
        end

    end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Create special experimental map
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    function Diag_special

        % Diagnostic map to show the image standard deviation in every
        % single voxel of the image relative to the mean of every voxel in
        % the image.

        % This is helpful to evaluate the temporal variability in an image.

        idat_std = std(single(idat),0,4);
        idat_map = single(imean)./idat_std;
        inparg.map_txt = 'Special Map';
        fsb_functional_map(idat_map,inparg,sandbox);

        disp('Special map shows the mean of every voxel in the image divided by ');
        disp('the standard deviation in every single voxel of the image');
        disp('This is helpful to evaluate the temporal stability of an image.');

    end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Create Sandbox trial difference map
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    function Diag_Sandmap

        % Diagnostic map to show the difference between trial and movement periods.
        % This is helpful to look for timing errors in the dataset
        if ~isempty(sandbox.intrial);
            findtrial = find(sandbox.intrial(:,1)==1);
            findnotrial = find(sandbox.intrial(:,1)==0);
            inparg.sandmap = mean(inparg.idat(:,:,:,findtrial),4)-mean(inparg.idat(:,:,:,findnotrial),4);
        end

        idat_map = inparg.sandmap;
        inparg.map_txt = 'Trial difference Map';
        fsb_functional_map(idat_map,inparg,sandbox);

        disp('Trial difference map shows the difference between trial and movement periods.');

    end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Create AB experimental map
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    function Diag_AB
        % Diagnostic map to show parts of the image whose time courses
        % correlate especially well with noise or artefacts in areas outside of the brain.
        % Shown is the correlation coefficient between the parts of
        % an image that are outside of the brain and the whole image. Areas
        % with a high correlation coefficient are most probably artefacts.
        % This is useful to determine the amount of ghosting in an image.

        idat_single = single(idat);

        slicey = int16(round(b/2));
        slicez = int16(round(c/2));

        roi_out1 = idat(:,slicey+round(slicey/1.5):b,...
            1:slicez+round(slicez/3),:);%outside brain area back

        roi_out2 = idat(:,1:round(slicey/4),...
            1:slicez+round(slicez/3),:);%outside brain area front

        roi_out = [roi_out1 roi_out2]; % combine both areas
        h = waitbar(0,'Calculating AB map...');
        for x = 1:d;
            waitbar((x/10)/d);
            roi_avg(x,1) = fsb_avg(roi_out(:,:,:,x));
        end

        idat_map = zeros(a,b,c);

        for y=1:b;
            waitbar(y/b);
            for x=1:a;
                for z =1:c;
                    R = corrcoef(squeeze(idat_single(x,y,z,:)),roi_avg);
                    idat_map(x,y,z) = R(2,1);
                end
            end
        end
        close(h);
        inparg.map_txt = 'Andrei Belitski Experimental Map';
        fsb_functional_map(idat_map,inparg,sandbox);

        disp('The AB map shows the correlation between the timecourse of a mean value')
        disp('of a predefined fixed area outside of the brain and the voxels in the image.')
        disp('Values will tend to get large if there is strong ghosting')
        disp ('or if the outside area includes part of the brain.')
        disp('This is helpful to evaluate ghosting artefacts in a scan.');

    end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Create advanced experimental map
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    function Diag_Ghost

        % Ghost map:
        % show parts of the image whose time courses
        % correlate especially well with noise or artefacts in areas outside of the brain.
        % Shown is the correlation coefficient between the parts of
        % an image that are outside of the brain and the whole image. Areas
        % with a high correlation coefficient are most probably artefacts.
        % This is useful to determine the amount of ghosting in an image.

        %~~~~~~~~~~~~~
        % This part is essentially the same as in the special map
        %~~~~~~~~~~~~~
        idat_single = single(idat);
        idat_std = std(idat_single,0,4);
        idat_map = idat_std./single(imean);

        %~~~~~~~~~~~~~
        % Find the means and standard deviations for the special map
        %~~~~~~~~~~~~~
        map_mean = mean(idat_map(:));
        map_std = std(idat_map(:));

        %~~~~~~~~~~~~~
        % Use the values to determine a threshold
        %~~~~~~~~~~~~~
        map_thr = map_mean+map_std;
        map_mask = ones(size(idat_map));
        map_mask(idat_map<map_thr) = NaN;
        map_mask = repmat(map_mask,[1 1 1 d]);
        roi_out = idat_single.*map_mask;

        h = waitbar(0,'Calculating Ghost map, please wait...');

        %~~~~~~~~~~~~~
        % Determine roi average
        %~~~~~~~~~~~~~
        roi_avg = zeros(d,1);
        for x = 1:d;
            waitbar((x/10)/d);
            single_roi = roi_out(:,:,:,x);
            roi_avg(x,1) = nanmean(single_roi(:));
        end

        %~~~~~~~~~~~~~
        % Calculate correlation between ROI and idat voxels
        %~~~~~~~~~~~~~
        for y=1:b;
            waitbar(0.1+(y/1.1)/b);
            for x=1:a;
                for z =1:c;
                    R = corrcoef(squeeze(idat_single(x,y,z,:)),roi_avg);
                    idat_map(x,y,z) = R(2,1);
                end
            end
        end
        close(h);

        inparg.map_txt = 'Ghost Map';
        fsb_functional_map(idat_map,inparg,sandbox);
        disp('The ghost map shows the correlation between the timecourse of a mean value')
        disp('of the area outside of the brain determined by a custom algorithm')
        disp('and every voxel in the image.')
        disp('Values will tend to get large in areas in which ghosting is likely to have an impact.')
        disp('This is helpful to evaluate ghosting artefacts in a scan.');

    end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Create SNR map
% This divides the mean of every voxel within the image by the
% Standard deviation of a predefined noise area outside of the brain
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    function Diag_SNR

        [a b c d] = size(idat);

        
%         inparg.noisetd = [size(idat,3)-3 size(idat,3)];
%         noise = single(idat(inparg.noise(1):inparg.noise(2),inparg.noise(3):inparg.noise(4),inparg.noisetd(1):inparg.noisetd(2),inparg.slice_n(4)));
        d = 1;

        if d>1 % for functional data
            inparg.noise2 = nanstd(single(inparg.noise),0,4);
        else % for single volumes
            inparg.noise2 = nanstd(inparg.noise(:));
        end

        noise_std = nanmean(inparg.noise(:));

        % if you want to have the whole image stack
        idat_map = single(imean)./noise_std;
        % if you want only the single image you have selected
        idat_map = single(idat(:,:,:,inparg.slice_n(4)))./noise_std;

        idat_map(isnan(idat_map))=0;
        idat_map(isinf(idat_map))= 0;
        inparg.map_txt = 'SNR Map';
        if ~isfield(inparg,'thr_high')
            inparg.thr_high = max(idat_map(:));
        end
        try
            fsb_functional_map(idat_map,inparg,sandbox);
        catch
            errmsg = lasterror;
            disp('Error in file:'); disp(errmsg.stack(1).file); disp('Line:');disp(errmsg.stack(1).line);
            disp('SNR calculation not possible, brain already extracted?')
        end

        disp('SNR map divides the mean of every voxel within the image by the ')
        disp('Standard deviation of a predefined noise area outside of the brain');

    end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Create Trial number map
% This determines the amount of observations (Trials) that you probably will
% need in every voxel of the brain to get statistical significance at the
% level of 0.05, assuming thatyou have a 1% signal change

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    function Diag_TN

        [a b c d] = size(idat);
        istd = std(single(idat),0,4);
        imean2 = single(imean);
        idat_map = zeros(a,b,c);
        trialnum = max(sandbox.intrial(:,3));
        triallength = round(length(sandbox.intrial)/trialnum);
        signal_change = 1.01;
        h = waitbar(0,'Calculating sample size map, please wait...');
        for x = 1:a;
            waitbar(x/a);
            for y = 1:b
                for z = 1:c;
                    idat_map(x,y,z) = sampsizepwr('z',[imean2(x,y,z) istd(x,y,z)],imean2(x,y,z)*signal_change);
                end
            end
        end
        close (h);
        idat_map(idat_map>3*mean(imean2(:)))=0;
        idat_map = idat_map/(triallength/2);
        inparg.map_txt = 'Sample Size Map';
        fsb_functional_map(idat_map,inparg,sandbox);

        disp('Trial number map determines the amount of observations (Trials) that ');
        disp('you probably will need in every voxel of the brain to get');
        disp('statistical significance at the level of 0.05, assuming that');
        disp('you have a 1% signal change');

    end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Create STD map
% Standard deviation map displays the Standard Deviation of every
% single voxel in the brain
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    function Diag_STD

        try
            idat_map = std(single(idat),0,4);
        catch
            for x=1:size(idat,1)
                idat_slice = (single(idat(x,:,:,:)));
                idat_map(x,:,:) = std(idat_slice,0,4);
            end;
        end

        inparg.map_txt = 'Standard Deviation Map';
        inparg.thr_high = max(idat_map(:));
        inparg.thr_low = min(idat_map(:));
        fsb_functional_map(idat_map,inparg,sandbox);

        disp('Standard deviation map displays the Standard Deviation of every')
        disp('single voxel in the brain')

    end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Create STD map normalized to image intensity
% Standard deviation map displays the Standard Deviation of every
% single voxel in the brain
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    function Diag_STD2

        idat_map = std(single(idat),0,4)./imean;
        inparg.map_txt = 'Standard Deviation Map normalized to signal intensity';
        fsb_functional_map(idat_map,inparg,sandbox);

        disp('Standard deviation map normalized displays the Standard Deviation of every')
        disp('single voxel in the brain normalized to its average intensity')

    end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Create CV map
% The Coefficient of Variance(CV) map displays the STD of every voxel over time
% divided by the mean of every voxel in the brain
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    function Diag_CV

        idat_map = std(single(idat),0,4);
        idat_map = idat_map./single(imean);
        inparg.map_txt = 'CV Map';
        fsb_functional_map(idat_map,inparg,sandbox);

        disp('Coefficient of Variance(CV) map displays the STD of every voxel over time')
        disp('divided by the mean of every voxel in the brain')

    end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Compute z_score map
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

% The timepoints of every single voxel in a trial are subtracted from the mean
% and the result is then divided by the standard deviation of every
% voxel

    function Diag_zscore

        idat_temp = permute(idat,[4 1 2 3]);
        idat_map = zscore(single(idat_temp));
        idat_map = permute(idat_map,[2 3 4 1]);
        idat_map = mean(idat_map,4);
        inparg.map_txt = 'z-score Map';
        fsb_functional_map(idat_map,inparg,sandbox);

        disp('Zscore map displays mean zscore in every voxel in the brain')
    end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Create volume histogram
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

% get a histogram of all voxel values in the selected volume

    function Diag_Hist

        idat_hist = double(idat(:,:,:,inparg.slice_n(4)));
        idat_hist = idat_hist(:);
        idat_hist(idat_hist==0) = NaN;
        figure;
        hist(idat_hist,100);

        disp('Histogram of voxel values')
    end

end
